library(plyr)
library(estimatr)
library(emmeans)
library(tidyverse)
library(magrittr)
library(broom)
library(specr)
library(metafor)

d1 <- "https://github.com/thomasjwood/full_fact/raw/main/ff_mt.rds" %>% 
  url %>% 
  gzcon  %>% 
  readRDS %>% 
  mutate(
    out_agree = out_agree %>% 
      mapvalues(
        c("Tend to disagree",
          "Tend to agree"),
        c("Disagree",
          "Agree")
      ) %>% 
      factor(
        c("Strongly disagree",
          "Disagree",
          "Neither agree nor disagree",
          "Agree",
          "Strongly agree")
      ),
    out_true = out_true %>% 
      factor(
        c("True",
          "Probably true",
          "Not sure",
          "Probably false",
          "False") %>% 
          rev
      ),
    agree_num = out_agree %>% 
      as.numeric,
    true_num = out_true %>% 
      as.numeric,
    cond = cond %>% 
      factor(
        c("cond_items",
          "cond_misinfo",
          "cond_corr")
      )
  ) %>% 
  filter(
    wave == "wave 1"
    ) %>% 
  mutate(
    gender = gender %>% 
      is_in(
        c("Female",
          "Male")
      ) %>% 
      ifelse(
        gender, NA
      ),
    com_num = agree_num %>% 
      add(
        true_num
      ) %>% 
      divide_by(2),
    country = country %>% 
      str_detect("saltwater|cooling") %>% 
      ifelse(
        "global", country
      ),
    employ2 = employ %>% 
      mapvalues(
        c("Retired",
          "Full-time parent, homemaker",
          "Self employed", 
          "Student/Pupil",
          "Employed full-time",
          "Employed part-time", 
          "Unemployed and not looking for a job/Long-term sick or disabled", 
          "Unemployed but looking for a job"),
        c(NA,
          "Unemployed",
          "Employed", 
          "Student",
          "Employed",
          "Employed", 
          "Unemployed", 
          "Unemployed")
      ) %>% 
      factor(
        c("Employed",
          "Unemployed",
          "Student")
      ),
    educ2 = educ %>% 
      mapvalues(
        
        c("NVQ4 / HNC / HND / Bachelor's degree or similar.",
          "Secondary", 
          "University",
          "High school completed",
          "Complete University", 
          
          "NVQ3/ SCE Higher Grade/ Advanced GNVQ/ GCE A/AS or similar.", 
          "University degree completed",
          "Technikon diploma/degree completed", 
          "NVQ5 or post-graduate diploma.",
          "GNVQ / GSVQ / GCSE/ SCE standard.",
          
          "Complete Secondary",
          "Incomplete University",
          "Artisans certificate, technical/secretarial qualification obtained", 
          "Complete Tertiary education (non-University)",
          "Post graduate degree completed", 
          
          "Incomplete Tertiary education (non-University)",
          "Complete Postgraduate / Master / Doctor's Degree", 
          "NVQ1, NVQ2",
          "Higher non-university",
          "Secondary school (age under 15 years old)", 
          
          "Some high school",
          "Incomplete Postgraduate / Master / Doctor's Degree", 
          "Incomplete Secondary",
          "Primary",
          "Complete Primay",
          
          "Incomplete Primary", 
          "No formal education",
          "Some primary school"),
        
        c("Bachelor's",
          "Secondary", 
          "Bachelor's",
          "Secondary",
          "Bachelor's", 
          
          "Bachelor's", 
          "Bachelor's",
          "Post secondary", 
          "Post secondary",
          "Secondary",
          
          "Secondary",
          "Post secondary",
          "Post secondary", 
          "Post secondary",
          "Master's or higher", 
          
          "Post secondary",
          "Master's or higher", 
          "Secondary",
          "Post secondary",
          "Secondary", 
          
          "Primary",
          "Master's or higher", 
          "Primary",
          "Primary",
          "Primary",
          
          "Primary", 
          "Primary",
          "Primary")
      ) %>% 
      mapvalues(
        c("Bachelor's",
          "Master's or higher",
          "Primary",
          "Secondary"),
        c("Bachelor's or higher",
          "Bachelor's or higher",
          "Secondary or less",
          "Secondary or less")
      ) %>% 
      factor(
        c("Secondary or less",
          "Post secondary",
          "Bachelor's or higher")
      ),
    ideo2 = ideo %>% 
      mapvalues(
        c("Left1", "1", "2", "3", "4",
          "5", "6",
          "7", "8", "9",  "10", "Right10", 
          
          "DK", "Don't know", "Prefer not to say"),
        c("Left", "Moderate", "Right", "DK") %>% 
          rep(
            c(5, 2, 5, 3)
          )
      ) %>% 
      factor(
        c("Left", "Moderate", "Right", "DK")
      ),
    age = 2020 %>% 
      subtract(
        birthyr
      ) %>% 
      cut(
        breaks = c(-Inf, 26, 40, 60, Inf),
        labels = c(
          "18-26", "27-40", "41-60", "41-60"
        )
      ) 
  )

t1 <- d1 %>% 
  group_by(country, issue) %>% 
  nest %>% 
  full_join(
    tibble(
      covar = c(
        "gender",
        "ideo2",
        "employ2",
        "educ2",
        "age"
        )
    ),
    by = character()
  )



t1$em <- map2(
  t1$data,
  t1$covar,
  function(i, j)
    
    lm_robust(
      "com_num ~ cond *" %>% 
        str_c(j) %>% 
        as.formula,
      data = i
      ) %>% 
      emmeans(
        "consec ~ cond | " %>%
          str_c(j) %>%
          as.formula
        )
  )

t1$covar %<>%
  mapvalues(
    c(
      "gender",
      "ideo2",
      "employ2",
      "educ2",
      "age"
    ),
    c("Gender", "Ideology", "Employment", "Education", "Age")
  )

t2 <- t1 %>% 
  ungroup %>% 
  select(
    country, issue, covar, em
  ) %>% 
  pmap_dfr(
    \(country, issue, covar, em){
    
      j <- em %>% 
        extract2("contrasts") %>% 
        tidy
      
      names(j)[[1]] <- "covar_x"
      
      j %>% 
        mutate(
          country, issue, covar
        )
      }
  ) %>% 
  mutate(
    contrast = contrast %>% 
      str_detect(" - cond_items") %>% 
      ifelse(
        "Misinformation effect",
        "Correction effect"
      ) %>% 
      factor(
        c("Misinformation effect",
          "Correction effect")
        )
    ) %>% 
  group_by(
    covar, 
    covar_x,
    contrast
  ) %>% 
  nest

t2$ma <- t2$data %>% 
  map(
    \(i)
    rma(
      yi = i$estimate, 
      sei = i$std.error
    )
  )


t3 <- t2 %>% 
  ungroup %>% 
  select(
    data, ma, covar, covar_x, contrast
    ) %>% 
  pmap_dfr(
    \(
      data, ma, covar, covar_x, contrast
      )

    ma %>% 
      tidy %>% 
      mutate(covar, covar_x, contrast)
    )


t4 <- t2 %>% 
  ungroup %>% 
  pmap_dfr(
    \(covar_x, contrast, covar, data, ma)
    
    data %>% 
      mutate(
        covar_x, contrast, covar
        )
    )
    

t4$covar_x %<>%
  factor(
    t4$covar_x %>% 
    unique %>% 
      rev
    )

t3$lab <- t3$estimate %>% 
  round(2) %>% 
  as.character %>% 
  str_replace("0\\.", ".")

t3$lab <- t3$estimate  %>%  
  is_less_than(0) %>% 
  ifelse(
    t3$lab %>% 
      str_pad(width = 4, side = "right", pad = "0"),
    t3$lab %>% 
      str_pad(width = 3, side = "right", pad = "0")
    )

t5 <- t4 %>% 
  group_by(covar, contrast) %>% 
  nest

t5$stat <- t5 %>% 
  use_series(data) %>% 
  map_dbl(
    function(i)
      aov(i$estimate ~ i$covar_x) %>% 
      tidy %>% 
      extract2("statistic") %>% 
      extract(1)
    )

t4$covar %<>%
  factor(
    t5 %>% 
      group_by(covar) %>% 
      summarize(mu = stat %>% mean) %>% 
      arrange(
        desc(mu)
        ) %>% 
      use_series(covar)
    )

t3$covar %<>%
  factor(t4$covar %>% levels)

t4 %>% 
  ggplot(
    aes(
      covar_x, estimate
    )
  ) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .5) +
  geom_dotplot(
    binaxis = "y",
    stackdir ='center',
    stackratio =.65, 
    dotsize = .6,
    shape = 21,
    fill = "grey50",
    ) +
  geom_linerange(
    aes(
      covar_x, 
      ymin = estimate %>% subtract(std.error %>% multiply_by(1.96)),
      ymax = estimate %>% add(std.error %>% multiply_by(1.96))
    ),
    color  ="white",
    size = 3,
    data = t3
  ) +
  geom_linerange(
    aes(
      covar_x, 
      ymin = estimate %>% subtract(std.error %>% multiply_by(1.96)),
      ymax = estimate %>% add(std.error %>% multiply_by(1.96))
      ),
    size = .8,
    data = t3
  ) +
  geom_point(
    aes(
      covar_x, 
      y = estimate
      ),
    shape = 21,
    fill = "white",
    size = 6.5,
    data = t3
  ) +
  geom_text(
    aes(
      covar_x, 
      y = estimate,
      label = lab
    ),
    size = 2.5,
    data = t3,
    family = "Roboto"
  ) +
  coord_flip() +
  facet_grid(
    covar ~ contrast,
    scales = "free",
    space = "free"
    ) +
  labs(
    y = "Difference on 5pt scale",
    x = ""
  ) +
  scale_x_discrete(
    breaks = t4$covar_x %>% 
      levels,
    labels = t4$covar_x %>% 
      levels %>% 
      str_replace_all(
        c("Bachelor's or higher" = "Bachelor's\nor higher" ,
          "Post secondary" = "Post\nsecondary", 
          "Secondary or less" = "Secondary\nor less")
      )
  ) +
  theme(
    strip.text.x = element_text(
      # size =  11, 
      angle = 0
    ),
    strip.text.y = element_text(
      # size =  11, 
      angle = 0
    ),
    legend.position = "none",
    axis.text.y = element_text(
      # angle = 45/2,
      hjust = .5,
      lineheight = .6
      # vjust = 1
    ),
    plot.margin = margin(.2, .2, .2 , 1, "cm") 
  )

